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A  nodal  triangle-based  spectral  element  (SE)  method  for  the  shallow 
water  equations  on  the  sphere  is  presented.  The  original  SE  method  uses 
quadrilateral  elements  and  high-order  nodal  Lagrange  polynomials,  con¬ 
structed  from  a  tensor-product  of  the  Legendre-Gauss-Lobatto  points.  In 
this  work  we  construct  the  high-order  Lagrange  polynomials  directly  on 
the  triangle  using  nodal  sets  obtained  from  the  electrostatics  principle  [17] 
and  Fekete  points  [29].  These  points  have  good  approximation  properties 
and  far  better  Lebesgue  constants  than  any  other  nodal  set  derived  for  the 
triangle.  By  employing  triangular  elements  as  the  basic  building-blocks  of 
the  SE  method  and  the  Cartesian  coordinate  form  of  the  equations,  we  can 
use  any  grid  imaginable  including  adaptive  unstructured  grids.  Results  for 
six  test  cases  are  presented  to  confirm  the  accuracy  and  stability  of  the 
method.  The  results  show  that  the  triangle-based  SE  method  yields  the 
expected  exponential  convergence  and  that  it  can  be  more  accurate  than 
the  quadrilateral-based  SE  method  even  while  using  30%  to  60%  fewer 
grid  points  especially  when  adaptive  grids  are  used  to  align  the  grid  with 
the  flow  direction.  However,  at  the  moment,  the  quadrilateral-based  SE 
method  is  twice  as  fast  as  the  triangle-based  SE  method  because  the  latter 
does  not  yield  a  diagonal  mass  matrix. 
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1.  INTRODUCTION 

The  receut  treud  towards  distributed-iueiuory  coiuputers  having  thousands  of 
commodity  processors  has  rekindled  interest  in  the  development  of  local  high-order 
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methods  for  the  simulation  of  geophysical  fluid  dynamics  applications.  The  most 
common  local  high-order  method  is  the  spectral  element  (SE)  method.  The  SE 
method  can  be  constructed  in  modal  (spectral)  or  nodal  (physical)  space.  In  ad¬ 
dition,  the  building-blocks  (or  element  shapes)  used  for  the  SE  method  have  been 
the  quadrilateral  or  the  triangle.  If  quadrilaterals  are  used  then  the  SE  method  is 
typically  employed  in  nodal  space  where  the  basis  functions  are  constructed  from  a 
tensor-product  of  the  one-dimensional  Legendre  cardinal  functions  [16].  However, 
if  triangles  are  used  then  the  SE  method  has  been  typically  employed  only  in  modal 
space  [26].  This  was  due  to  the  lack  of  a  good  set  of  nodal  points  for  the  triangle. 

Using  the  electrostatics  principle,  Hesthaven  [17]  obtained  a  set  of  nodal  points 
on  the  triangle  with  good  approximation  properties  for  all  polynomials  of  order 

<  11.  Eor  11  <  <  15  we  use  the  Eekete  points  [29]  which  are  only  currently 

available  for  orders  N  <  20.  Using  this  nodal  set  Warburton  et  al.  [31]  showed 
exponential  convergence  for  the  incompressible  Navier-Stokes  equations.  The  suc¬ 
cess  of  these  results  has  inspired  us  to  seek  similarly  successful  applications  of  this 
nodal  triangular  set  for  the  solution  of  the  shallow  water  equations  on  the  sphere. 

The  shallow  water  equations  are  a  set  of  first  order  nonlinear  hyperbolic  equations 
which  contain  all  of  the  horizontal  operators  found  in  the  primitive  atmospheric 
equations  used  in  numerical  weather  prediction  (NWP)  and  climate  models.  Thus 
to  design  a  good  atmospheric  model  requires  a  good  shallow  water  model.  The 
construction  of  fast,  accurate,  and  flexible  atmospheric  models  is  the  ultimate  goal 
of  our  research.  In  this  quest,  we  have  successfully  developed  an  exponentially  con¬ 
vergent  global  atmospheric  model  using  the  nodal  quadrilateral-based  SE  method 
[13,  14].  While  the  accuracy  and  performance  of  this  model  have  been  shown  to 
exceed  those  of  operational  spectral  transform  models,  developing  adaptive  grids 
for  quadrilateral  elements  may  prove  too  cumbersome  to  pursue.  The  existence 
of  numerous  adaptive  triangular  mesh  generation  packages  [1,  8]  motivates  us  to 
explore  nodal  triangle-based  SE  methods. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2  describes  the  gov¬ 
erning  equations  of  motion  used  to  test  our  numerical  method.  In  Sec.  3  we  describe 
the  discretization  of  the  governing  equations.  This  includes  the  spatial  discretiza¬ 
tion  by  the  triangle-based  SE  method  and  the  time-integrator.  In  Sec.  4  we  describe 
a  few  of  the  many  possible  triangular  tessellations  of  the  sphere.  Einally,  in  Sec.  5 
we  present  convergence  rates  of  the  nodal  triangle-based  SE  method  and  compare 
it  with  the  quadrilateral-based  SE  method.  This  then  leads  to  some  conclusions 
about  the  feasibility  of  this  approach  for  constructing  future  atmospheric  models 
and  a  discussion  on  the  direction  of  future  work. 


2.  SHALLOW  WATER  EQUATIONS 

The  shallow  water  equations  are  a  system  of  first  order  nonlinear  hyperbolic 
equations  which  govern  the  motion  of  an  inviscid  incompressible  fluid  in  a  shallow 
depth.  The  predominant  feature  of  this  type  of  fluid  is  that  the  characteristic  length 
of  the  fluid  is  far  greater  than  its  depth  which  is  analogous  to  the  motion  of  air 
in  the  atmosphere  and  water  in  the  oceans.  Eor  this  reason  these  equations  are 
typically  used  as  a  first  step  toward  the  construction  of  either  NWP,  climate,  or 
ocean  models. 
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The  spherical  shallow  water  equations  in  Cartesian  advective  form  are 

^  +  V  •  {cf>u)  =  0  (1) 

— +u-Vu  = -f  {x  X  u) -V  +  -  lix,  (2) 

where  the  nabla  operator  is  defined  as  V  =  {dx,dy,dz)^ ,  cj)  is  the  geopotential 
height  {(j)  =  gh  where  g  is  the  gravitational  constant  and  h  is  the  vertical  height 
of  the  fluid),  is  the  surface  topography  (e.g.,  mountains),  u  =  {u,v,w)'^  is 
the  Cartesian  wind  velocity  vector,  /  =  2^  is  the  Coriolis  parameter  and  (w,a) 
represent  the  rotation  of  the  earth  and  its  radius,  respectively. 

The  term  fix,  where  x  =  {x,y,z)'^  is  the  position  vector  of  the  grid  points,  is  a 
fictitious  force  introduced  to  constrain  the  fluid  particles  to  remain  on  the  surface 
of  the  sphere.  By  switching  from  spherical  (2D)  to  Cartesian  (3D)  coordinates  we 
have  allowed  the  fluid  particles  an  additional  degree  of  freedom  which  will  manifest 
itself  in  the  fluid  particles  flying  off  the  surface  of  the  sphere.  In  order  to  prevent 
this  from  happening  we  introduce  the  Lagrange  multiplier  /i. 

The  shallow  water  equations  in  Cartesian  form  have  received  significant  attention 
recently  (see  [28,  7,  9,  10,  15,  11,  20]).  It  should  be  mentioned  that  the  Cartesian 
form  of  the  equations  introduces  no  approximations  whatsoever;  the  equations  are 
completely  equivalent  to  the  equations  in  spherical  coordinates  as  shown  by  Swarz- 
trauber  et  al.  [28].  The  reason  for  using  the  Cartesian  form  of  the  equations  is  that 
the  pole  singularity  associated  with  spherical  coordinates  is  avoided  and  because 
this  form,  in  conjunction  with  the  SE  mapping  described  in  Sec.  3.1.2,  allows  for 
any  curved  surface  to  be  discretized  by  this  approach.  For  example,  we  could  easily 
change  the  shape  of  the  spherical  domain  to  any  warped  spheroidal.  This  will  allow 
for  a  more  realistic  representation  of  the  earth  and  other  planets  which  are  not 
perfect  spheres. 


3.  DISCRETIZATION 

In  this  section  we  describe  the  discretization  of  the  shallow  water  equations. 
In  Sec.  3.1,  we  describe  the  spatial  discretization  by  the  spectral  element  method 
including:  the  choice  of  basis  functions,  integration,  and  the  Alter  used  to  maintain 
stability.  In  Sec.  3.2  we  discuss  the  explicit  time-integrator,  the  time-filter  used  to 
control  the  computational  modes,  and  the  Lagrange  multiplier  required  to  constrain 
the  fluid  particles  onto  the  sphere. 

3.1.  Spatial  Discretization  by  the  Triangle-based  Spectral  Element 

Method 

3.1.1.  Background  and  Motivation 

To  explain  the  need  for  the  electrostatics  and  Fekete  points  on  the  triangle  a  few 
words  regarding  interpolation  theory  are  required.  Let  Pjv  denote  the  space  of  all 
polynomials  of  degree  <  N  and  be  a  collection  of  points  in  Pjy.  In  addition,  let 
q  be  an  arbitrary  function  that  we  wish  to  interpolate  and  In  an  unique  function 
in  Pjv  such  that  7jv(9(^i))  =  q{^i)-  Using  the  usual  definition  of  the  Loo  norm 

II  q  11=  maxjeQ  |  q{^)  \ 
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II  In  ||—  niax||g||=i  |  In{q)  \ 

we  can  now  quantify  the  error  of  approximating  q  by  In{q)  within  the  triangle  fi. 
Assuming  the  existence  of  a  function  r  G  Pjv  which  best  approximates  q  and  noting 
that  we  can  write  r  =  we  can  compute  the  interpolation  error  of  q  as 

II  <1  -  lN{q)  11=11  q-r  +  lN{r)  -  lN{q)  II  • 

Using  the  Cauchy-Schwarz  inequality  yields 

II  q  -  lN{q)  ||<  (1+  II  In  id  \\q-r\\  (3) 

where  ||  In  ||  denotes  the  Lebesgue  constant  which  must  be  minimized  in  order  to 
avoid  the  well-known  Runge  effect  [26]  typical  of  monomial  expansions  of  the  type 
I  0  <  i  +  y  <  N .  This  is  the  expansion  typically  associated  with  the  standard 
finite  element  method  where  the  points  (^,7?)  are  equi-spaced  within  the  trian¬ 
gle.  Since  it  is  not  known  how  to  compute  Lebesgue  points  (i.e.,  points  which 
minimize  the  Lebesgue  constant)  we  use  points  computed  via  an  electrostatics 
analogy  or  points  which  maximize  the  determinant  of  the  Vandermonde  matrix 
(Fekete  points)  which  indirectly  approximate  them.  On  the  1-simplex  (the  line 
—  1  <  ^  <  1)  the  points  satisfying  the  electrostatics  and  Fekete  principles  are  in 
fact  the  Legendre-Gauss-Lobatto  (LGL)  points.  However,  on  the  2-simplex  (the  tri¬ 
angle  — 1<^,  7?<1;  ^  +  7?<0)  determining  which  are  the  optimal  points  remains 
an  open  question.  Nonetheless,  attempts  have  been  made  to  construct  nodal  sets 
which  yield  LGL  points  along  the  edges  of  the  triangle.  Both  the  electrostatics  and 
Fekete  points  satisfy  this  condition  and  the  fact  that  they  have  good  Lebesgue  con¬ 
stants  make  them  suitable  choices  for  constructing  triangle-based  spectral  element 
basis  functions.  The  Fekete  and  electrostatics  points  coincide  with  the  standard 
equi-spaced  points  for  N  <2.  However,  the  Lebesgue  constant  for  the  Fekete  (VA^) 
and  electrostatics  points  {N  <  9)  increases  proportionally  with  N  while  for  the 
standard  equi-spaced  points  the  Lebesgue  constant  increases  exponentially  with  N. 

The  term  ||  q  —  r  ||  in  Eq.  (3)  is  minimized  if  the  polynomial  space  Pjv  is  able 
to  approximate  the  function  q.  We  use  the  Proriol-Koornwinder-Dubiner  (PKD) 
polynomials  [24,  19,  5]  as  they  are  analogous  to  the  ID  Legendre  polynomials  with 
the  properties  that  they  are  stably  computable  and  are  a  conveniently  orthogonal 
basis  defined  on  the  triangle.  The  construction  of  the  cardinal  functions  based  on 
the  PKD  polynomials  is  the  topic  of  the  next  section.  In  this  section,  we  have  only 
discussed  the  interpolation  properties  of  certain  nodal  sets  and  bases,  however,  in 
order  to  construct  discrete  operators  requires  not  only  good  interpolation  properties 
but  also  good  integration  properties. 

On  the  1-simplex  the  LGL  points  happen  to  be  both  optimal  interpolation  and 
integration  points  which  results  in  a  diagonal  mass  matrix  since  the  interpolation 
and  integration  points  are  co-located;  having  a  diagonal  mass  matrix  is  important 
for  achieving  efficiency.  On  the  2-simplex  such  points  have  not  yet  been  found  and 
thus  far  one  must  be  content  to  choose  either  good  interpolation  or  integration  but 
not  both.  In  this  paper  we  choose  good  interpolation  and  then  use  exact  numerical 
integration  formulas  (cubature)  which,  while  quite  accurate,  does  not  result  in  a 
diagonal  or  lumpable  mass  matrix,  since  the  interpolation  and  integration  points 
are  not  co-located. 
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It  should  be  noted  that  the  Fekete  points  have  been  used  as  integration  points 
in  previous  works  [18,  23,  30]  for  linear  problems  which  then  results  in  a  diagonal 
mass  matrix;  however,  the  Fekete  points  only  integrate  up  to  order  N  polynomials 
and  so  for  nonlinear  terms  (quadratics  terms  such  as  2N)  this  order  of  integration 
is  insufficient  to  achieve  exponential  convergence.  Our  results  using  this  approach 
proved  to  be  quite  dismal.  For  this  reason  we  choose  to  use  2N  cubature  rules  with 
the  ill-effect  of  having  to  contend  with  a  global  sparse  mass  matrix.  We  should  also 
mention  that  our  initial  results  with  the  diagonal  mass  matrix  nodal  set  proposed  by 
Mulder  [22]  have  been  quite  encouraging;  however,  these  points  only  exist  for  <  5 
and  they  are  not  nearly  as  stable  as  the  approach  we  present  in  this  manuscript. 


3.1.2.  Basis  Functions 

To  define  the  local  operators  which  shall  be  used  to  construct  the  global  approx¬ 
imation  of  the  solution  we  begin  by  decomposing  the  spherical  domain  Q,  into 
non-overlapping  triangular  elements  ffe  such  that 

iVe 

0  =  U  Oe. 

e=l 

To  perform  differentiation  and  integration  operations,  we  introduce  the  nonsin¬ 
gular  mapping  a;  =  ^(4)  which  defines  a  transformation  from  the  physical  Carte¬ 
sian  coordinate  system  x  =  {x,y,z)'^  to  the  local  reference  coordinate  system 
4  Vj  C)  such  that  (^,  t})  lies  on  the  spherical  surface  tiled  by  the  triangular 
elements  defined  by  fig  =  {(Cj  vX)j  “1  <  ^  ^  ^  <  Oj  C  =  I}- 

Let  us  now  represent  the  local  element-wise  solution  q  by  an  A^th  order  polyno¬ 
mial  in  4  as 

Mjv 

k=l 

where  represents  Mjv  =  |(A^-|-l)(A^-|-2)  grid  points  and  are  the  associated 
multivariate  Lagrange  polynomials.  For  the  grid  points  we  choose  the  nodal 

set  derived  from  the  electrostatics  principle  [17]  for  A^  <  11  and  the  Fekete  points 
for  11<N  <  15. 

We  construct  the  Lagrange  polynomials  on  the  reference  triangle,  Lk{^,  rj),  which 
are  implicitly  defined  by  their  cardinal  nature,  by  reference  to  an  easily  constructed 
orthonormal  PKD  polynomial  basis  [24,  19,  5].  This  basis  is  defined  as 


(2i  +  l)(i+i  +  l) 

p0,0 


(^) 


(.v) 


(4) 


where  Pn’^iO  represents  the  nth  order  Jacobi  polynomial  in  the  interval  —1  < 
^  <  1,  k  =  i  +  j{N  -h  1)  -h  1,  and  the  indices  vary  as  0  <  +  j  <  N,  and 

fc  =  1, . . . ,  Mjv. 

We  next  seek  an  explicit  formula  for  the  Lagrange  basis  by  representing  them  in 
terms  of  the  reference  basis,  i.e. 


Mjv 

Li  itv)  =  '^Aikipk 

k=l 


(5) 
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where  the  indices  are  now  defined  as  i,j,k=  1, . . . ,  Mjy.  We  then  use  the  cardinal 
property  of  the  Lagrange  polynomials 


k=l 

where  6  is  the  Kronecker  delta  function,  to  determine  that 

Ak  =  ■  (6) 

We  next  recognize  that 


Vjk  —  4’k 

is  a  generalized  Vandermonde  matrix  and  using  Eqs.  (5),  (6),  and  (7)  we  construct 
the  Lagrange  polynomials  as  follows 


=  (8) 

k=l 

We  note  that,  for  example,  the  Vandermonde  matrix  generated  using  this  basis  and 
the  15th  order  Fekete  nodes  only  has  a  condition  number  of  approximately  116. 
This  good  conditioning  allows  us  to  use  this  approach  in  contrast  to  conventional 
wisdom  regarding  poor  conditioning  of  the  standard  monomial  based  Vandermonde. 

3.1.3.  Filtering  the  High-Frequency  Waves 

Recall  that  the  local  element-wise  Vandermonde  matrix  can  be  written  as 


Vij  —  Ipj  ■ 

Thus  the  function  values  at  each  grid  point  i  inside  the  element  fig  can  be  defined, 
using  a  modal  (spectral)  expansion,  as  follows 


qi  =  (9) 

1=1 

where  q  represents  the  expansion  coefficients  in  modal  space  of  the  function  q.  In 
matrix  form  we  can  write  Eq.  (9)  as 


q  =  V  q  (10) 

which  we  call  the  nodal  triangle-based  SE  transform  because  it  allows  us  to  trans¬ 
form  from  nodal  to  modal  space.  Inverting  Eq.  (10)  yields 

q  =  (V“^)^  q 

which  yields  the  amplitudes  in  the  modal  representation  (amplitude- frequency  space) 
We  can  then  filter  these  amplitudes  in  any  manner  but  here,  based  on  past  experi¬ 
ence  [10, 13],  we  choose  the  Boyd-Vandeven  transfer  function  [2]  which  we  denote  by 
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A.  Applying  the  filter  to  the  amplitudes  and  then  transforming  to  nodal  (physical) 
space  is  achieved  in  the  following  matrix-vector  multiply  operation 

F  C  y, 

q  =  F  q 


where 

F  =  VA  (V-i)^  (11) 

is  the  Mjv  x  Mjv  filter  matrix.  This  filter  matrix  is  applied  every  10  time-steps  at 
only  20  %  strength  (see  [10]). 

3.1.4-  Integration 

In  order  to  complete  the  discussion  of  the  local  element-wise  operations  required 
to  construct  discrete  spectral  element  operators  we  must  lastly  describe  the  inte¬ 
gration  procedure  required  by  the  weak  formulation  of  all  Galerkin  methods.  For 
any  two  functions  /  and  g  the  integration  X  proceeds  as  follows 

.  Mq 

m,9]  =  /  f{x)g{x)dx  =  I  J(4J  I  f{ii)g{ii) 

where  Mq  is  a  function  of  Q  which  represents  the  order  of  the  cubature  approxi¬ 
mation.  For  Wi  and  4*  we  use  the  high-order  cubature  rules  for  the  triangle  given 
in  [27,  3,  21,  4]  of  order  2N.  This  order  is  chosen  for  two  reasons:  first,  it  is  a 
good  compromise  between  accuracy  and  efficiency;  and  second,  it  allows  for  a  fair 
comparison  of  our  new  triangle-based  SE  method  with  the  quadrilateral-based  SE 
method  which  typically  uses  order  2N  —  1  quadrature  rules. 


3.1.5.  Local  Element-wise  Operators 

To  simplify  the  description  of  the  numerical  algorithm,  let  us  define  the  following 
local  element  operators:  let 


Mfj  =  /  Li{x)Lj{x)dx,  =  /  Li{x)Lj{x)VLk{x)dx, 

-D-J  =  /  Li{x)VLj{x)dx,  Li{x)Lj{x)Lk{x)Li{x)dx 


(12) 

(13) 


represent  the  mass,  advection,  differentiation,  and  Coriolis  matrices  where  i,j,  k,  I  = 
l,...,Mjv.  Note  that  A®  =  (A^,  A^,  A^)  and  =  (T>®,T>®,T>|)  are  vectors  of 
matrices  corresponding  to  the  three  spatial  directions. 


3.1.6.  Satisfying  the  Equations  Globally 

To  satisfy  the  equations  globally  requires  summing  the  local  element  matrices, 
Eqs.  (12)  and  (13),  to  form  their  global  representation.  This  summation  procedure 
is  known  as  the  global  assembly  or  direct  stiffness  summation.  Let  us  represent  this 
global  assembly  procedure  by  the  summation  operator 

iVe 

A 

e=l 
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with  the  mapping  (i,e)  — >  (I)  where  i  =  are  the  local  element  grid 

points,  e  =  are  the  spectral  elements  covering  the  spherical  shell,  and 

I  =  1, ...  ,Np  are  the  global  grid  points.  Applying  the  global  assembly  operator  to 
the  local  element  matrices  results  in  the  following  global  matrices: 

iVe  iVe  iVe  iV^ 

M=/\M®,  A=  f\A^,  D=  f\D^,  C=f\C^. 

e=l  e=l  e=l  e=l 

With  these  operators  defined  and  by  denoting  the  global  grid  vector  for  the  grid 
points  as  x,  the  geopotential  as  (j),  and  the  wind  velocity  as  u  we  can  now  state  the 
variational  form  of  the  problem  as:  find  {(j),  G  'i  L  £  such  that 

+  {A(l>f  u  =  -  (!>  (14) 

M ^  +  (Au)^  u  =  —Cf  {x  y.  u)  —  D  {(j)  +  (j)^)  —  M (iix)  (15) 

CJT 

where  JJ^(fl)  is  the  space  of  all  functions  with  functions  and  first  derivatives  be¬ 
longing  to  I/^(fl)  -  the  space  of  all  functions  that  are  square  integrable  over  fl. 
For  (j)  and  u  we  choose  the  polynomial  space  Pjv-Pjv  without  violating  the  inf-sup 
condition. 


3.2.  Time-Integrator  and  the  Lagrange  Multiplier 

Discretizing  the  equations  in  time  by  the  leapfrog  method  yields 

-  2At  ((A^)^  u  +  A^u(^  "  (16) 

-  2 At  {{Auf  u  +  Cf{x  x  it))"  (17) 


Since  we  need  to  ensure  that  the  velocity  field  remains  tangential  to  the  sphere, 
we  require 

a;  •  It  =  0. 

Let  us  first  write  Eq.  (17)  as 

it"+i  =  -  2Atiix,  (18) 

where 

=  Mit"-^  -  2At  (^{Auf  It  -F  Cf{x  X  It))"  -  2At  {D{(j)  +  ^®))"  . 
Taking  the  scalar  product  of  Eq.  (18)  with  x,  and  rearranging  yields 
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Using  Eq.  (19)  we  note  that  we  can  project  any  3D  vector  q  onto  the  surface  of  the 
sphere  by  applying  the  following  matrix  operator 

PQ  =  q-  ^ix  q)x  (20) 

where  the  projection  matrix  P  is  given  by 

—xy  —xz  \ 

P  =  —  i  -xy  -  y^  -yz  J  .  (21) 

\  —xz  —yz  c?  —  z^  ) 

This  allows  the  equations  in  constrained  form  to  be  written  as  follows 

-  2At  u  +  A^ucf^  "  (22) 

=  PB^.  (23) 

Since  M  is  a  global  sparse  symmetric  matrix  we  use  the  conjugate  gradient  method 
with  point  Jacobi  preconditioning  to  invert  it.  Static  condensation  is  used  to  sep¬ 
arate  the  boundary  and  interior  points  of  the  elements  in  order  to  reduce  the  com¬ 
putational  cost  of  inverting  M. 

3.2.1.  Time-Filter 

To  complete  the  discussion  of  the  time- integrator,  we  must  describe  the  Asselin 
filter  used  for  the  leapfrog  method  which  otherwise  would  admit  a  nonphysical 
mode  to  propagate  in  the  opposite  direction  from  the  physical  mode.  This  filter  is 
applied  as  follows 

q"  =  q"  +  y(q"+i-2r+g"^^)  (24) 

where  the  tilde  represents  the  unfiltered  solution  and  7  =  0.01  is  the  filter  strength. 
While  it  is  true  that  for  7  7^  0  the  leapfrog  method  becomes  first  order,  we  have 
used  this  time-integrator  because  it  is  typically  used  in  atmospheric  and  shal¬ 
low  water  models  which  will  then  facilitate  comparisons  with  other  methods.  For 
more  suitable  time- integrators  for  the  shallow  water  and  atmospheric  equations  the 
reader  is  referred  to  the  recent  papers  [11,  12,  14]  where  semi-Lagrangian,  operator- 
integration-factor  splitting,  and  backward  difference  methods  are  discussed. 

4.  GRID  GENERATION  ON  THE  SPHERE 

The  choice  of  which  triangulation  to  use  for  the  sphere  is  not  obvious.  Commonly, 
grids  are  chosen  which  simplify  the  construction  of  the  discrete  operators.  For 
example,  latitude-longitude  grids  are  used  with  spectral  transform  methods  because 
these  are  the  only  grids  that  can  be  used  with  this  method.  The  hexahedral  grid 
(i.e.,  the  cubed-sphere)  has  been  used  with  finite  difference,  finite  volume,  and 
spectral  element  methods  because  each  of  the  six  faces  of  the  cube  map  onto  a 
simple  Cartesian  geometry  that  allows  for  the  simple  and  rapid  construction  of  the 
discrete  operators.  Picking  one  grid  and  constructing  the  discrete  operators  on  a 
specific  grid  geometry  simplifies  matters  but  it  also  dictates  the  algorithm  thereby 
losing  any  hope  of  using  other  types  of  grids  and  adaptive  solution  strategies. 
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In  our  case,  the  spectral  element  method  is  constructed  in  a  very  general  way 
such  that  the  model  reads  in  any  grid  geometry  and  then  constructs  the  discrete 
operators  directly  on  the  grid.  This  allows  the  use  of  any  grid  and  offers  the  free¬ 
dom  to  choose  the  best  possible  grid  for  specific  applications.  For  the  purposes 
of  validating  the  triangle-based  SE  method  we  shall  use  a  disjointed  set  of  trian¬ 
gles  formed  by  the  subdivision  of  the  triangular  faces  of  an  icosahedron;  however, 
it  should  be  understood  that  any  triangular  grid  can  be  used.  We  shall  compare 
the  triangle-based  SE  method  with  the  quadrilateral-based  SE  method  described 
in  [9,  11,  13].  With  the  quadrilateral-based  SE  method  we  used  hexahedral  and 
icosahedral  grids.  The  choice  of  hexahedral  grids  is  a  natural  one  because  it  repre¬ 
sents  a  structured  grid  that  has  become  quite  popular  with  many  newly  proposed 
grid  point  methods.  The  quadrilateral-based  icosahedral  grid  is  chosen  because  it 
is  an  unstructured  grid  and  represents  a  good  comparison  for  the  triangle-based 
icosahedral  grid.  Einally,  to  illustrate  the  geometric  flexibility  of  the  triangle-based 
SE  method  we  describe  a  number  of  triangular  unstructured  grids  formed  by  the 
Delaunay  triangulation  of  various  point  sets. 


4.1.  Triangle-based  Icosahedral  Grids 

To  construct  icosahedral  grids  we  consider  the  initial  icosahedron  and  subdivide 
each  of  the  initial  triangles  by  a  Lagrange  polynomial  of  order  nj.  Prior  to  mapping 
these  elements  onto  the  sphere  it  is  convenient  to  map  the  triangles  onto  a  gnomonic 
space.  The  most  unbiased  mapping  is  obtained  by  mapping  about  the  centroid  of 
the  triangles. 

Let  (Ac,¥?c)  be  the  centroid  of  the  triangle  we  wish  to  map  where  A  represents 
the  zonal  (east-west)  and  y?  the  meridional  (north-south)  directions.  The  gnomonic 
mapping  is  then  given  by 

^  ^  _ Qcosy?sin(A- Ac) _ 

sin  tpc  sin  (p  +  cos  Pc  cos  p  cos(A  —  Ac)  ’ 

^  a  [cos  Pc  sin  p  —  sin  p^  cos  p  cos(A  —  Ac)] 
sin  Pc  sin  p  cos  Pc  cos  p  cos(A  —  Ac) 

where  X  G  [—1,  -1-1]^  in  the  equi-distant  gnomonic  space  G.  To  simplify  matters  a 
bit,  we  first  apply  a  rotation  whereby  Eq.  (25)  becomes 


X  =  atanAi{,  E  =  a  tan  sec  A^j  , 


in  the  new  coordinate  system  with  the  centroid  {Xc,Pc)  located  at  (0,0).  The 
rotation  mapping  is  given  as 


fR 


arctan 


cos¥?sin(A  —  Ac) 


sinyjc  sinyj  -I-  cos  Pc  cospcos{X  —  Ac)  J 
arcsin  [cos  Pc  sin  p  —  sin  pc  cos  p  cos(A  —  Ac)] . 


This  approach  enables  the  construction  of  a  triangle-based  icosahedral  grid  with 
the  following  properties 


Np  =  W{niNf+2 
Nc  =  20{nif  , 


(26) 

(27) 
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where  Np  and  denote  the  number  of  points  and  elements  comprising  the  trian¬ 
gular  grid,  and  nj  controls  the  number  of  triangular  elements  while  N  denotes  the 
order  of  the  polynomial  inside  each  element.  Figure  1  provide  examples  of  grids  for 


a)  b)  c) 


FIG.  1.  Triangle- based  icosahedral  grid  for  nj  =  2  and  a)  A  =  4,  b)  A  =  8,  and  c)  A  =  12. 

nj  =  2  and  N  =  4,8  and  16.  All  the  grids  illustrated  are  viewed  from  the  North 
Pole  where  the  thick  lines  denote  the  elements  and  the  thin  lines  me  the  high-order 
grid  points. 

Since  we  use  the  icosahedral  grid  to  compare  the  triangle-based  SE  method  with 
the  standard  quadrilateral-based  SE  method  on  icosahedral  and  hexahedral  grids 
it  is  important  to  compare  grids  with  comparable  total  grid  points.  Erom  [13]  the 
number  of  points  for  the  quad-based  icosahedral  and  hexahedral  grids  are  given  as 

Np  =  60{nf  Nf  +  2  (28) 

and 

Np  =  6{n%  Nf  +  2  (29) 

where  nf  and  are  parameters  which  controls  the  number  of  elements  in  the 
grids.  Equating  Eqs.  (26),  (28),  and  (29)  shows  that  the  triangle-based  icosahedral 
grid  has  the  same  number  of  grid  points  as  the  quadrilateral-based  grids  for 

m  =  V6nf  =  n%  (30) 

which  we  approximate  by  rij  ~  2nj  ~  Thus  rij  =  3,  nj  =  2  yield 

approximately  the  same  number  of  grid  points. 

It  should  be  mentioned  that  these  grid  constructions  do  not  actually  yield  the 
same  number  of  grid  points  and  elements.  In  fact  for  this  set  of  grid  constructions 
the  triangle-based  icosahedral  grid  has  80  triangular  elements  while  the  quadrilateral- 
based  icosahedral  and  hexahedral  grids  have  60  and  54,  respectively.  Although  for 
a  given  polynomial  order  {N)  they  all  span  the  same  polynomial  space  within  each 
element  they  differ  greatly  in  the  total  number  of  grid  points  and  in  their  distri¬ 
bution.  In  fact,  the  triangle-based  icosahedral  grid  has  33%  fewer  points  than  the 
quad-based  icosahedral  grid  and  25%  fewer  points  than  the  quad-based  hexahe¬ 
dral  grid.  In  addition,  the  triangle-based  grid  appears  to  be  far  more  isotropic 
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than  the  quadrilateral  grids.  In  other  words,  the  triangle-based  icosahedral  grid 
has  no  biasing  in  its  orientation  and  thereby  treats  all  directions  equally.  Unlike 
the  quadrilateral-based  icosahedral  grid,  the  quadrilateral-based  hexahedral  grid  is 
structured  and  isotropic  but  only  along  two  directions  in  each  of  the  six  faces  of  the 
hexahedron.  Therefore,  the  hexahedral  grid  will  do  quite  well  for  flows  which  are 
aligned  with  these  two  directions.  We  shall  see  that  some  of  the  test  cases  benefit 
from  this  type  of  grid  orientation  while  others  do  not. 

4.2.  Triangulations  Based  on  the  Platonic  Solids 

In  the  previous  section  we  described  the  construction  of  triangular  grids  using  the 
icosahedron.  However,  there  is  nothing  special  about  this  Platonic  solid;  we  could 
have  used  any  of  the  five  Platonic  solids.  For  example,  one  could  use  the  tetrahedron 
(4  triangular  faces),  hexahedron  (6  quadrilateral  faces),  octahedron  (8  triangular 
faces),  dodecahedron  (12  pentagonal  faces),  or  icosahedron  (20  triangular  faces). 
However,  the  faces  of  the  hexahedron  and  dodecahedron  need  to  be  subdivided 
into  triangular  elements  (24  for  the  hexahedron  and  60  for  the  dodecahedron). 
Using  the  Platonic  solids  one  can  create  a  grid  with  the  following  properties 

Nv  =  ^{nNf+2 
Ne  =  NtU^ 

where  Nt  are  the  number  of  initial  triangles  comprising  the  Platonic  solid  and 
n  is  the  subdivision  of  each  of  the  initial  triangles.  Nt  =  4,  8,  20,  24,  and  60  for 
the  tetrahedron,  octahedron,  icosahedron,  hexahedron,  and  dodecahedron,  respec¬ 
tively.  In  Fig.  2  we  show  equivalent  resolution  triangulations  using  the  octahedron, 
icosahedron,  and  hexahedron  viewed  from  the  North  Pole.  From  Fig.  2  it  is  evident 


a) 


b) 


c) 


FIG.  2.  Triangulations  of  the  sphere  using  the  a)  octahedron  {Ne  =  64),  b)  icosahedron 
[Ne  =  80),  and  c)  hexahedron  {Ne  =  96)  viewed  from  the  North  Pole. 


that  the  octahedral  and  icosahedral  grids  are  more  uniform  than  the  hexahedral 
grid;  this  is  due  to  the  fact  that  the  initial  triangles  of  the  octahedron  and  icosahe¬ 
dron  are  equilateral  whereas  those  for  the  hexahedron  are  isosceles.  However,  there 
are  hardly  any  differences  between  these  sets  of  grids  but  due  to  the  popularity  of 
icosahedral  grids  we  shall  use  this  set  in  Sec.  5  but  it  should  be  understood  that  all 
these  grids  yield  similar  convergence  rates. 
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4.3.  Triangulations  Based  on  Latitude-Longitude 

Because  triangles  are  the  2-simplex  this  means  one  can  construct  an  unique  tri¬ 
angulation  for  random  point  sets  distributed  on  the  surface  of  the  sphere.  This 
then  facilitates  the  construction  of  triangulations  on  the  sphere  because  as  long  as 
one  can  create  point  sets  they  can  be  triangulated  quite  naturally  using  Delaunay 
triangulation  methods.  For  the  latitude- longitude  grids  we  use  the  STRIPACK 
Delaunay  triangulation  software  [25]. 

In  this  section  we  discuss  a  few  of  the  many  possible  triangular  grids  that  can  be 
constructed  on  the  sphere.  We  limit  the  grids  to  fixed  adaptive  grids;  that  is,  grids 
that  are  constructed  at  the  beginning  of  the  time-integration  and  remain  fixed  for 
all  time.  In  the  future  we  shall  address  the  dynamically  adaptive  approach. 

Latitude-longitude  grids  are  perhaps  the  simplest  tessellations  for  the  sphere.  In 
fact,  they  are  popular  for  NWP  and  climate  models  because  these  are  the  only 
grids  that  spectral  transform  methods  can  use.  Figure  3a  shows  a  regular  latitude- 
longitude  grid  with  {Nion,Niat)  =  (40,20).  On  each  latitude  ring  (East- West  di¬ 
rection)  this  grid  has  the  same  number  of  points  {Nion  =  40)  and  the  grid  spacing 
becomes  smaller  as  we  approach  the  poles.  This  situation  is  known  as  the  pole 
problem  because  there  are  too  many  redundant  points  near  the  poles  as  is  illus¬ 
trated  by  the  lines  of  constant  longitude  (vertical)  as  they  curve  toward  the  poles. 
In  addition,  the  proximity  of  adjacent  points  near  the  poles  severely  restricts  the 
maximum  time-step  that  can  be  used.  One  way  around  this  dilemma  is  to  use  thin 
or  reduced  grids.  As  we  approach  the  poles,  the  number  of  points  on  each  latitude 
ring  is  decreased  in  order  to  maintain  a  constant  grid  spacing  throughout.  This 
type  of  grid  is  depicted  in  Fig.  3b  where  lines  of  constant  longitude  no  longer  curve 
toward  the  poles.  Both  the  regular  and  thin  latitude-longitude  grids  can  be  used 


FIG.  3.  Triangulations  of  the  sphere  using  the  following  types  of  latitude- longitude  grids: 
a)  regular  {Ne  =  1600),  b)  thin  {Ne  =  1044),  and  c)  adaptive  {Ne  =  1844)  viewed  from  (0,0). 

with  the  spectral  transform  method  (without  the  triangular  elements).  However, 
if  the  discretization  method  is  based  on  triangles,  then  we  can  go  even  further  and 
refine  the  region  of  interest  and  use  a  coarse  grid  everywhere  else.  Figure  3c  shows 
an  adaptive  grid  where  the  region  of  interest  is  along  the  Equator.  We  shall  use 
this  grid  to  show  the  advantages  that  adaptive  triangular  grids  may  offer. 
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5.  NUMERICAL  EXPERIMENTS 

For  the  numerical  experiments,  we  use  the  normalized  I/2  error  norm 


UWl. 


/^(^exact  -  4>?d^ 
Sq  ^exact 


to  judge  the  accuracy  of  the  SE  methods.  To  compute  the  Courant  number  the 
elements  are  decomposed  into  their  high-order  (HO)  grid  points  and  these  grid 
points  form  a  fine  grid  which  we  refer  to  as  the  HO  cells.  The  velocities  and  grid 
spacings  are  then  defined  at  the  centers  of  these  cells.  Using  these  definitions  the 
Courant  number  is  then  defined  as 


Courant  number  =  max 


VeG  [l,...,7Ve] 


where 


J  U  for  case  1 

\  -h  for  cases  2, 3, 4, 5,  and  6 


where  C  is  the  characteristic  speed,  U  =  y'u  ■  u,  and  As  =  -|-  Ay^  +  Az^  is 

the  grid  spacing.  For  all  the  results  presented  the  Courant  number  is  taken  to  be 

0.5. 

Six  test  cases  are  used  to  judge  the  performance  of  the  triangle-based  SE  method. 
Cases  1,  2,  3,  5,  and  6  correspond  to  the  Williamson  et  al.  standard  test  case  suite 
[32].  Case  4  was  recently  proposed  by  Galewsky  et  al.  for  testing  shallow  water 
models  [6].  This  case  presents  a  more  challenging  test  than  those  in  the  Williamson 
et  al.  test  suite  because  if  the  resolution  is  not  sufficiently  high  then  the  numerics 
will  not  be  able  to  sustain  the  steady  zonal  jet  with  steep  vorticity  gradient.  If  the 
method  cannot  sustain  the  jet  then  the  accuracy  declines  rapidly.  Case  1  involves 
the  geopotential  equation  (passive  advection)  only  whereas  the  remainder  of  the 
test  cases  concern  the  full  nonlinear  shallow  water  equations.  In  addition,  cases  1, 
2,  3,  and  4  have  analytic  solutions  whereas  cases  5  and  6  do  not  and  are  only  used 
to  determine  the  accuracy  of  the  triangle-based  SE  method  qualitatively.  These 
last  two  test  cases  have  been  run  by  a  vast  community  and  the  results  are  well- 
documented  for  comparison. 


5.1.  Description  of  the  Test  Cases 

5.1.1.  Case  1:  Passive  Advection  of  a  Cosine  Wave 

Case  1  concerns  the  solid  body  rotation  of  a  cosine  wave.  The  velocity  field 
remains  unchanged  throughout  the  computation.  Williamson  et  al.  [32]  recommend 
that  the  error  be  computed  after  12  days  of  integration  which  corresponds  to  one 
complete  revolution  of  the  cosine  wave. 

5.1.2.  Case  2:  Steady-State  Nonlinear  Zonal  Ceostrophic  Flow 

This  case  is  a  steady-state  solution  to  the  nonlinear  shallow  water  equations.  The 
equations  are  geostrophically  balanced  and  remain  so  for  the  duration  of  the  inte¬ 
gration  where  the  velocity  field  remains  constant  throughout  the  computation.  The 
geopotential  height  (p  undergoes  a  solid  body  rotation  but  since  the  initial  height 
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field  is  given  as  a  constant  band  in  the  zonal  direction  and  the  fiow  field  is  purely 
zonal,  then  the  solution  remains  unchanged  throughout  the  time-integration.  The 
velocity  field  is  the  same  as  that  used  in  case  1.  Williamson  et  al.  [32]  recommend 
that  the  error  be  computed  after  5  days  of  integration. 

5.1.3.  Case  3:  Steady-State  Nonlinear  Zonal  Geostrophic  Flow  with  Compact 
Support 

This  case  is  another  steady-state  solution  to  the  nonlinear  shallow  water  equa¬ 
tions  where  the  equations  remain  geostrophically  balanced  for  the  duration  of  the 
integration.  The  initial  velocity  field  is  zero  everywhere  except  in  a  very  small 
isolated  region.  This  isolated  region,  or  jet,  encapsulates  the  fiow  and  confines  the 
geopotential  height  field  to  remain  within  a  localized  circular  region.  The  results 
are  reported  for  a  5-day  integration  as  suggested  in  [32]. 

5.1.4.  Case  4-  Galewsky  et  al.  Zonal  Dynamics 

This  test  case  consists  of  a  zonal  jet  and  an  unperturbed  balanced  initial  geopo¬ 
tential  height  field.  The  balanced  initial  field  should  be  maintained  indefinitely  but 
Galewsky  et  al.  [6]  suggest  running  the  case  for  5  days.  This  is  a  rather  stringent 
test  of  shallow  water  models  because  if  the  accuracy  and/or  the  resolution  is  not 
sufficiently  high  then  the  model  will  not  be  able  to  sustain  the  balanced  initial  field 
and  the  error  will  increase  quite  rapidly,  unlike  cases  1,  2,  and  3  which  are  much 
more  forgiving.  In  addition,  because  the  jet  is  zonally  positioned,  then  any  grid  that 
is  not  aligned  with  the  zonal  direction  will  have  much  more  difficulty  maintaining 
the  jet. 

5.1.5.  Case  5:  Zonal  Flow  over  an  Isolated  Mountain 

This  case  uses  the  same  initial  conditions  as  case  2  with  the  addition  of  a  conical 
mountain  at  (A,yj)  =  (180,30).  Due  to  the  zonal  fiow  impinging  on  the  mountain, 
various  wave  structures  form  and  propagate  throughout  the  sphere.  This  test  is 
run  for  15  days  as  suggested  in  [32]. 

5.1.6.  Case  6:  Rossby-Haurwitz  Wave 

Although  Rossby-Haurwitz  waves  are  not  analytic  solutions  to  the  shallow  water 
equations,  they  have  become  a  de  facto  test  case.  In  a  non-divergent  barotropic 
model,  the  initial  geopotential  height  field  undergoes  a  solid  body  rotation  in  a 
counterclockwise  direction  when  viewed  from  the  North  Pole.  Although  this  case 
does  not  have  an  analytic  solution,  it  is  well-known  that  the  initial  wave  structure 
of  the  Rossby-Haurwitz  wave  should  remain  intact  for  the  duration  of  the  time- 
integration. 


5.2.  Results  on  the  Test  Case  Suite 

5.2.1.  Convergence  Rate  of  the  Triangle-based  SE  Method 

Before  analyzing  the  behavior  of  the  triangle-based  SE  method  on  the  standard 
test  case  suite  or  delving  into  comparisons  between  the  triangle-based  SE  method 
with  the  quadrilateral-based  method  it  is  important  to  quantify  the  convergence 
rate  of  the  method.  Eor  this  study  we  use  case  2  and  show  convergence  rates  in 
Eig.  4  for  A^  <  14.  The  order  of  convergence  shown  in  Eig.  4  is  computed  as  an 
average  convergence  rate  computed  over  all  the  grid  refinements  where  at  each  grid 
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FIG.  4.  Case  2.  The  normalized  </>  L2  error  as  a  function  of  icosahedral  grid  refinement, 
ni,  after  5  days  for  the  polynomial  orders  N  =  1,2,4,6,8,10,12,  and  14  with  their  associated 
convergence  rates. 


refinement  the  convergence  rate  is  defined  as 

log  [error„^ + 1  /error„^  ] 
i^ate  =  — ri — f — ^ 

log  [ni/{ni  +  1)J 

This  figure  shows  that  the  expected  spectral  accuracy  is  achieved  for  all  values  of 
N,  that  is, 

error  oc 

which  states  that  the  error  decreases  exponentially  with  increasing  N. 

In  [18,  30,  23]  the  Fekete  points  are  used  as  both  interpolation  and  integration 
points.  The  advantage  of  this  strategy  is  that  the  triangle-based  SE  method  yields 
a  diagonal  mass  matrix.  An  obvious  question  to  ask  is  why  have  we  chosen  a 
different  approach,  that  is,  to  use  the  Fekete/Electrostatics  points  for  interpolation 
and  the  high-order  cubature  rules  from  [27,  3,  21,  4]  for  integration  which  then  no 
longer  yields  a  diagonal  mass  matrix.  To  answer  this  question  we  run  case  2  for 
the  icosahedral  grid  m  =  2  and  =  8  for  30  days.  The  results  of  this  simulation 
are  shown  in  Eig.  5  where  it  is  clear  that  the  2N  cubature  rule  which  we  use 
throughout  this  paper  is  far  more  accurate  and  stable  for  long  time-integrations 
than  the  N  cubature  rule.  Thus  while  the  N  cubature  rule  may  be  appropriate  for 
steady-state  problems  or  linear  time-dependent  problems  it  is  not  appropriate  for 
nonlinear  problems  with  long  time-integrations  such  as  those  typical  in  geophysical 
fluid  dynamics  applications. 

5.2.2.  Case  1:  Passive  Advection  of  a  Cosine  Wave 

Figure  6  shows  that  the  SE  method  on  triangles  and  quadrilaterals  converge 
algebraically  regardless  of  the  structure  of  the  grid.  Exponential  convergence  is 
not  expected  for  this  case  because  the  derivative  at  the  base  of  the  cosine  hill  is 
non-smooth.  Note  that  the  triangle-based  and  quadrilateral-based  SE  converge  at 
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FIG.  5.  Case  2.  The  normalized  </>  L2  error  as  a  function  of  time  in  days  for  the  2N  cubature 
rule  (solid  line)  which  does  not  yield  a  diagonal  mass  matrix  and  the  N  cubature  rule  (dashed 
line)  which  does  yield  a  diagonal  mass  matrix  for  nj  =  2  and  N  =  8. 


approximately  the  same  rate  with  the  quadrilateral-based  method  being  slightly 
better  than  the  triangle-based  method. 


FIG.  6.  Case  1.  The  normalized  (p  L2  error  as  a  function  of  polynomial  order,  N,  after 
12  days  using  80,  60,  and  54  elements  for  the  icosahedral  triangle-based  (Ico  Tri),  icosahedral 
quadrilateral- based  (Ico  Quad),  and  hexahedral  quadrilateral- based  (Hex  Quad)  spectral  elements. 


5.2.3.  Case  2:  Steady-State  Nonlinear  Zonal  Geostrophic  Flow 
Figure  7  illustrates  that  the  triangle-based  and  quadrilateral-based  SE  methods 
yield  exponential  convergence.  The  quadrilateral-based  hexahedral  grid  yields  the 
best  results  with  the  triangle-based  icosahedral  grid  giving  better  results  than  the 
quadrilateral-based  icosahedral  grid.  The  hexahedral  grid  will  be  very  difhcult  to 
beat  for  this  test  case  because  it  is  aligned  with  the  direction  of  the  flow.  For  this 
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reason  it  is  expected  that  a  latitude-longitude  grid  would  also  do  extremely  well  for 
this  test  case. 


FIG.  7.  Case  2.  The  normalized  </>  L2  error  as  a  function  of  polynomial  order,  N,  after  5  days 
using  80,  60,  and  54  elements  for  the  icosahedral  triangle- based  (Ico  Tri),  icosahedral  quadrilateral- 
based  (Ico  Quad),  and  hexahedral  quadrilateral-based  (Hex  Quad)  spectral  elements. 


5.2.4-  Case  3:  Steady-State  Nonlinear  Zonal  Geostrophic  Flow  with  Compact 
Support 

Figure  8  shows  that  the  triangle-based  and  quadrilateral-based  SE  methods, 
again,  yield  exponential  convergence;  however,  the  triangle-based  SE  method  is 
slightly  better  than  the  quadrilateral-based  SE  methods.  Unlike  case  2,  this  case 
represents  a  more  localized  flow  problem.  The  circular  region  where  the  geopoten¬ 
tial  is  confined  is  not  aligned  with  any  of  the  grids  and  so  the  isotropy  of  the  grid 
becomes  an  important  factor.  This  explains  why  the  triangle-based  SE  method  on 
the  icosahedral  grid  yields  better  accuracy  than  the  quadrilateral-based  methods. 

5.2.5.  Case  4-  Galewsky  et  al.  Zonal  Dynamics 

Eigure  9  illustrates  the  convergence  rates  of  the  triangle-based  and  quadrilateral- 
based  SE  methods  as  a  function  of  total  number  of  grid  points,  Np,  with  8th  order 
polynomials.  This  figure  shows  that  triangle-based  and  quadrilateral-based  SE 
methods  yield  high-order  convergence  and  that  the  triangle-based  SE  method  is 
superior  to  the  quadrilateral-based  SE  methods.  Once  again,  the  absence  of  grid 
alignment  with  the  flow  has  hindered  the  quadrilateral-based  SE  methods  from 
resolving  the  physics  of  the  problem  as  well  as  the  more  isotropic  triangle-based 
icosahedral  grid. 

Eigures  10  and  11  show  snapshots  of  the  grid,  geopotential  (^),  and  zonal  velocity 
{us)  after  a  5-day  integration  for  the  triangle-based  and  quadrilateral-based  SE 
methods  with  various  grid  resolutions  using  8th  order  polynomials.  The  view  of  the 
figures  is  from  the  North  Pole.  In  Eig.  10  it  is  observed  that  the  resolutions  nj  =  4 
{Np  =  10242)  and  below  are  not  capable  of  maintaining  the  jet.  The  formation  of 
wave  number  5  structures  begin  to  emerge  due  to  the  lack  of  symmetry  in  the  grid 
near  the  initial  points  of  the  icosahedron  (which  form  a  pentagon  around  the  pole) . 
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FIG.  8.  Case  3.  The  normalized  </>  L2  error  as  a  function  of  polynomial  order,  N,  after  5  days 
using  80,  60,  and  54  elements  for  the  icosahedral  triangle- based  (Ico  Tri),  icosahedral  quadrilateral- 
based  (Ico  Quad),  and  hexahedral  quadrilateral-based  (Hex  Quad)  spectral  elements. 


FIG.  9.  Case  4.  The  normalized  (p  L2  error  as  a  function  of  grid  points,  Np  after  5 
days  using  N  =  8  for  the  icosahedral  triangle-based  (Ico  Tri),  icosahedral  quadrilateral-based  (Ico 
Quad),  and  hexahedral  quadrilateral-based  (Hex  Quad)  spectral  elements. 


For  m  =  5  {Np  =  16002)  the  balanced  initial  state  is  maintained  where  the  jet  is 
visibly  intact.  Similarly,  Fig.  11  shows  that  the  quadrilateral-based  SE  method  on 
the  hexahedral  grid  cannot  maintain  the  jet  for  resolutions  Np  <  38402.  Thus  the 
triangle-based  SE  method  on  icosahedral  grids  can  maintain  the  jet  with  58%  fewer 
points  than  the  quadrilateral-based  SE  method  on  hexahedral  grids. 

5.2.6.  Case  5:  Zonal  Flow  over  an  Isolated  Mountain 

Eigure  12  shows  snapshots  of  the  fields  after  a  15-day  integration  for  the  triangle- 
based  SE  method  with  nj  =  6  and  N  =  8  and  the  quad-based  SE  with  nn  =  8 
and  N  =  8.  The  view  of  the  figures  is  from  (A,¥?)  =  (180,0).  Note  that  the 
number  of  grid  points  between  the  two  methods  is  approximately  equal.  In  fact. 
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b) 

FIG.  10.  Case  4.  Contours  of  the  grid  (left),  4>  (center),  and  Us  (right)  for  the  triangle-based 
SE  method  after  5  days  for  A  =  8  and  a)  nj  =  4  {Np  =  10242)  and  b)  nj  =  5  {Np  =  16002) 
viewed  from  the  North  Pole. 


b) 

FIG.  11.  Case  4.  Contours  of  the  grid  (left),  4>  (center),  and  Us  (right)  for  the  quadrilateral- 
based  SE  method  after  5  days  for  A  =  8  and  a)  =  9  (Np  =  31106)  and  b)  =  10 
(Np  =  38402)  viewed  from  the  North  Pole. 


the  results  using  either  quadrilaterals  or  triangles  are  virtually  indistinguishable. 
All  grid  resolutions  beyond  this  one  yield  the  correct  wave  structures  for  both 
methods  with  little  or  no  differences  obtained  with  increased  resolution. 


5.2.7.  Case  6:  Rossby-Haurwitz  Wave 
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b) 

FIG.  12.  Case  5.  Contours  of  </>  (left),  Us  (center),  and  Vs  (right)  after  15  days  for  a)  the 
triangle-based  SE  method  with  nj  =  6  and  N  =  8  {Np  =  23042)  and  b)  the  quadrilateral-based 
SE  method  with  uh  =  8  and  N  =  8  {Np  =  24578)  viewed  from  (180,0)  which  is  in  the  vicinity  of 
the  mountain. 


Figure  13  shows  snapshots  of  the  fields  after  a  14-day  integration  for  the  triangle- 
based  SE  method  with  nj  =  6  and  N  =  8  and  the  quad-based  SE  with  uh  =  8  and 
N  =  8.  The  view  of  the  figures  is  from  the  North  Pole.  Once  again,  the  results 
of  the  triangle  and  quadrilateral  based  SE  methods  are  identical.  Eurthermore,  all 
grid  resolutions  beyond  this  one  yield  the  correct  wave  structures  for  both  methods 
with  minor  differences  obtained  with  increased  resolution. 

5.2.8.  The  Advantages  of  using  Adaptive  Grids 

To  show  the  advantages  that  adaptive  grids  may  offer  we  first  revisit  case  1.  Eor 
this  test  case  the  quadrilateral-based  hexahedral  grid  yielded  a  superior  convergence 
rate  to  the  triangle-based  SE  method  on  icosahedral  grids.  It  was  conjectured  that 
the  alignment  of  the  hexahedral  grid  with  the  direction  of  the  flow  accounted  for 
this  difference  in  convergence  rates.  Therefore,  let  us  now  compute  the  convergence 
rate  for  the  triangle-based  SE  method  on  an  adaptive  grid  aligned  with  the  flow 
direction.  Eigure  14  shows  the  grid  (with  N  =  8)  and  the  contours  of  (f)  after  one 
revolution.  Eigure  15  shows  that  the  the  convergence  rate  for  the  triangle-based  SE 
method  on  the  adaptive  grid  is  now  superior  to  that  of  the  quadrilateral-based  SE 
method  on  both  the  hexahedral  and  icosahedral  grids.  It  should  be  mentioned  that 
all  three  grids  have  approximately  the  same  number  of  grid  points.  Constructing 
this  type  of  grid  with  the  quadrilateral-based  SE  method  is  not  possible. 

As  a  second  example  of  the  advantages  of  adaptive  grids  we  revisit  case  4.  Recall 
that  with  the  isotropic  icosahedral  grid  it  required  16,000  grid  points  in  order  to 
maintain  the  jet  for  the  duration  of  the  5-day  integration.  In  Eig.  16  we  show  the 
grid,  (j),  and  u  contours  for  an  adaptive  grid  with  only  11,000  grid  points.  More 
importantly,  the  polynomial  order  is  =  6.  Thus  with  less  grid  points  and  with  a 
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b) 

FIG.  13.  Case  6.  Contours  of  4>  (left),  Us  (center),  and  Vs  (right)  after  14  days  for  a)  the 
triangle-based  SE  method  with  nj  =  6  and  N  =  8  (Np  =  23042)  and  b)  the  quadrilateral-based 
SE  method  with  njj  =  8  and  N  =  8  {Np  =  24578)  viewed  from  the  North  Pole. 


FIG.  14.  Case  1.  The  grid  (left)  and  contours  of  4>  (right)  for  the  triangle-based  SE  method 
for  N  =  8. 


lower  polynomial  order  the  adaptive  grid  solution  is  able  to  capture  and  maintain 
the  jet. 

As  a  final  example  of  adaptive  grids  we  revisit  case  5.  Recall  that  with  the 
isotropic  icosahedral  grid  it  required  23,000  grid  points  in  order  to  achieve  a  con¬ 
verged  solution  for  the  duration  of  the  15-day  integration.  In  Fig.  17  we  show 
the  grid,  (j),  u,  and  v  contours  for  an  adaptive  grid  with  only  18,000  grid  points. 
More  importantly,  the  polynomial  order  is  =  4.  Thus  with  less  grid  points  and 
with  a  lower  polynomial  order  the  adaptive  grid  solution  is  once  again  able  to  ob¬ 
tain  the  proper  wave  structure.  The  adaptive  grid  is  constructed  in  the  Northern 
Hemisphere  thereby  handling  the  wave  structure  quite  well  in  this  region.  Note, 
however,  that  due  to  the  coarser  grid  in  the  Southern  Hemisphere  some  of  the  wave 
structures  are  under-resolved. 
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FIG.  15.  Case  1.  The  normalized  </>  L2  error  as  a  function  of  polynomial  order,  N,  after 
12  days  using  92,  60,  and  54  elements  for  the  adaptive  latitude-longitude  triangle-based  (Adapt 
Tri),  icosahedral  quadrilateral-based  (Ico  Quad),  and  hexahedral  quadrilateral-based  (Hex  Quad) 
spectral  elements. 


FIG.  16.  Case  4.  Contours  of  the  grid  (left),  (p  (center),  and  Us  (right)  for  the  triangle-based 
SE  method  after  5  days  for  N  =  6  using  an  adaptive  grid  with  Ng  =  614  {Np  =  11054)  viewed 
from  the  North  Pole. 


FIG.  17.  Case  5.  Contours  of  the  grid  (left),  (p  (left-center),  Us  (right-center),  and  Vs  (right) 
for  the  triangle-based  SE  method  after  15  days  for  N  =  4  using  an  adaptive  grid  with  Ng  =  2308 
{Np  =  18466)  viewed  from  (180,0). 


6.  CONCLUSIONS 

The  newly  proposed  triangle-based  SE  method  exhibited  its  expected  exponential 
convergence.  Furthermore,  the  triangle-based  SE  method  was  shown  to  be  superior 
to  the  quadrilateral-based  methods  for  three  of  the  four  cases  having  analytic  so¬ 
lutions.  Cases  5  and  6  have  been  typically  used  to  represent  challenging  cases  but 
the  triangle-based  SE  method  was  capable  of  representing  the  underlying  physics 
of  these  problems  with  a  modest  grid  resolution.  However,  case  4  showed  that  for 
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crude  resolutions,  the  wave  structure  collapses  therefore  representing  a  good  addi¬ 
tion  to  the  Williamson  et  al.  shallow  water  test  case  suite.  For  cases  1  and  2,  the 
quadrilateral-based  SE  method  on  hexahedral  grids  was  shown  to  be  superior  to 
the  triangle-based  SE  method  on  icosahedral  grids.  These  cases,  however,  represent 
flows  which  are  aligned  with  the  hexahedral  grid  and  for  this  reason  it  is  difficult 
to  compete  with  the  quadrilateral-based  SE  method.  However,  with  the  assistance 
of  an  adaptive  grid,  the  triangle-based  SE  method  proved  to  be  superior.  Cases  3 
and  4  have  strong  local  characteristics  whereby  the  geopotential  height  is  confined 
within  a  small  circular  region  encapsulated  by  a  jet.  Eor  these  cases,  the  hexahedral 
grid  is  no  longer  aligned  with  the  predominant  features  of  the  flow.  It  is  conjectured 
that  the  isotropy  of  the  triangle-based  SE  method  on  the  icosahedral  grid  is  what 
allows  the  triangle-based  SE  method  to  achieve  a  better  convergence  rate  than  the 
quadrilateral-based  SE  methods.  The  success  of  the  triangle-based  SE  method  for 
these  two  test  cases  are  encouraging  because  it  implies  that  while  the  method  can 
accurately  handle  flows  with  strong  global  features  it  is  much  better  suited  for  han¬ 
dling  flows  with  significant  local  characteristics.  It  is  therefore  a  good  candidate 
for  use  with  adaptive  grids  as  was  shown  for  cases  1,  4,  and  5. 

It  should  be  mentioned  that  we  have  not  shown  performance  comparisons  be¬ 
tween  the  triangle-based  and  quadrilateral-based  methods.  At  the  present  the 
quadrilateral-based  SE  method  is  at  a  far  more  mature  stage  than  its  triangle- 
based  counterpart.  Unlike  the  quadrilateral-based  SE  method,  the  triangle-based 
SE  method  does  not  give  rise  to  a  diagonal  mass  matrix,  M,  and  thereby  requires 
inverting  a  sparse  global  matrix  at  each  time  step.  Even  though  the  triangle-based 
SE  method  requires  30%  to  60%  fewer  grid  points  than  its  quadrilateral-based  coun¬ 
terpart  to  achieve  the  same  accuracy,  the  inversion  of  a  global  matrix  significantly 
reduces  the  efficiency  of  the  triangle-based  method.  In  this  paper  we  have  only  pre¬ 
sented  the  triangle-based  SE  method  in  its  continuous  form,  in  the  future  we  shall 
explore  this  method  in  its  discontinuous  form  which  obviates  the  need  to  invert 
a  global  matrix.  We  will  also  consider  more  efficient  domain  decomposition  tech¬ 
niques  as  well  as  co-located  interpolation  and  integration  points  which  will  improve 
the  apparent  inefficiency  in  the  current  approach.  The  exponential  accuracy  and 
geometric  flexibility  of  the  triangle-based  SE  method  merits  further  investigations 
into  its  use  in  all  realms  of  geophysical  fluid  dynamics. 
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